#include "parafluidsmoothingkernel.h"
#include "parafluidsphforce.h"
#include "parafluidsphneighbor.h"

using namespace ParaFluid;

ForceUpdater::~ForceUpdater()
{
	/* does nothing */
}

ForceUpdater::ForceUpdater()
	: AtomVoxelizer::ReadWriteVisitor()
	, _SPH(0)
{
	/* does nothing */
}

void ForceUpdater::updateForces(SPH & sph)
{
	_SPH = &sph;
	_SPH->atomVoxelizer.acceptVisitor(*this);
}

void ForceUpdater::updateForcesParallel(SPH & sph, int indexFirst, int indexLast)
{
	_SPH = &sph;
	_SPH->atomVoxelizer.acceptVisitor(*this, indexFirst, indexLast);
}

void ForceUpdater::visit(Atom & atom)
{
	atom.currentForce = Vector(0, 0, 0);
	
	class ForceAccumulator
		: public NeighborhoodReadOnlyVisitor
	{
		
	private:
		
		SPH * _SPH;
		
		Atom * _atom;
		
	public:
		
		~ForceAccumulator()
		{
			/* does nothing */
		}
		
		ForceAccumulator(SPH & sph, Atom & atom)
			: NeighborhoodReadOnlyVisitor()
			, _SPH(&sph)
			, _atom(&atom)
		{
			/* does nothing */
		}
		
		void visitNeighbor(const Atom & neighbor)
		{
			if(neighbor.currentDensity != 0)
			{
				// Compute Pressure
				Vector distance(
					 _atom->currentPosition.x - neighbor.currentPosition.x,
					 _atom->currentPosition.y - neighbor.currentPosition.y,
					 _atom->currentPosition.z - neighbor.currentPosition.z
				);
				float pressureAtom = _atom->substance.pressureFactor * (_atom->currentDensity - _atom->substance.restDensity);
				float pressureNeighbor = neighbor.substance.pressureFactor * (neighbor.currentDensity - neighbor.substance.restDensity);
				float pressureForceFactor = (neighbor.substance.mass * (pressureAtom + pressureNeighbor)) / neighbor.currentDensity;
				Vector pressureForce = SmoothingKernel::gradientSpiky(_SPH->smoothingDistance(), distance) * pressureForceFactor;
				
				// Accumulate Pressure
				_atom->currentForce += pressureForce;
				
				// Compute Viscosity
				Vector velocityDifference(
					((neighbor.currentPosition.x - neighbor.previousPosition.x) - (_atom->currentPosition.x - _atom->previousPosition.x)) / _SPH->lastTimeStep(),
					((neighbor.currentPosition.y - neighbor.previousPosition.y) - (_atom->currentPosition.y - _atom->previousPosition.y)) / _SPH->lastTimeStep(),
					((neighbor.currentPosition.z - neighbor.previousPosition.z) - (_atom->currentPosition.z - _atom->previousPosition.z)) / _SPH->lastTimeStep()
				);
				float viscosityForceFactor = (_atom->substance.viscosityFactor * neighbor.substance.mass * SmoothingKernel::laplacianViscosity(_SPH->smoothingDistance(), distance))  /  neighbor.currentDensity;
				Vector viscosityForce = velocityDifference * viscosityForceFactor;
				
				// Accumulate Viscosity
				_atom->currentForce += viscosityForce;
			}
		}
		
	} accumulator(*_SPH, atom);
	
	accumulator.visitNeighborhood(*_SPH, atom);
}



















